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i-Q ' Abstract 

We present an original method for reconstructing a three-dimensional object having two spatial 
dimensions and one spectral dimension from data provided by the infrared slit spectrograph on board the 
Spitzer Space Telescope. During acquisition, the light flux is deformed by a complex process comprising 
four main elements (the telescope aperture, the slit, the diffraction grating and optical distortion) before 
■ it reaches the two-dimensional sensor. 

^ ' The originality of this work lies in the physical modelling, in integral form, of this process of data 

|H ' formation in continuous variables. The inversion is also approached with continuous variables in a semi- 

, parametric format decomposing the object into a family of Gaussian functions. The estimate is built in 

a deterministic regularization framework as the minimizer of a quadratic criterion. 
^ ' These specificities give our method the power to over-resolve. Its performance is illustrated using real 

^.f-^ I and simulated data. We also present a study of the resolution showing a 1.5-fold improvement relative 



to conventional methods. 



, Index Terms 
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inverse problems, bayesian estipmation, over-resolved imaging, spectral imaging, irregular sampled, 
interpolation, IRS Spitzer. 

I. Introduction 

Since the end of the 1970's, infra-red to millimetric observations of the sky from space have brought 
about a revolution in practically all fields of astrophysics. It has become possible to observe distant 
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galaxies, perform detailed physicochemical studies of interstellar matter. Observations in the far infra-red 
are now possible thanks to new types of sensors (Ge:Ga Si:As semiconductors and bolometer arrays). 
The properties of these new sensors encouraged the astrophysicists of the Institut d'Astrophysique 
Spatiale (IAS) to work with researchers at the Laboratoire des Signaux et Systemes (L2S) in order 
to develop suitable processing methods. The spectral imaging work presented here was carried out in 
the framework of this cooperative effort. The aim is to reconstruct an over-resolved object having two 
spatial dimensions (a, and one spectral dimension A. Data provided by the Infrared Spectrograph 
(IRS) 11] on board the american Spitzer Space Telescope launched in 2003 are used to illustrate our 
work. Several sets of two-dimensional data are delivered by the Spitzer Science Center (SSC), each 
set being the result of an acquisition for a given satellite pointing direction. The data were acquired 
using a slit spectrograph, the operation of which is described in detail in part [III This instrument is 
located in the focal plane of the telescope. When the telescope is pointed towards a region of the sky, 
the spectrograph slit selects a direction of space a. The photon flux is then dispersed perpendicularly to 
the slit direction with a diffraction grating. The measurement is made using a two-dimensional sensor. 
A signal containing one spatial dimension a and the spectral dimension A is thus obtained. The second 
spatial dimension /3 is obtained by scanning the sky (modifying telescope pointing). This scanning has 
two notable characteristics: 

• it is irregular, because the telescope control is not perfect; 

• it is, however, measured with sub-pixel accuracy (eighth of a pixel). 

In addition, for a given pointing direction, : the telescope optics, the slit width and the sensor integration 
limit the spatial resolution while the grating, the slit and the sensor integration limit the spectral resolution. 
The specificity of systems of this type is that the width of impulse response depends on the wavelength. 
A phenomenon of spectral also aliasing appears for the shortest wavelengths. Finally, the scanning results 
in irregular sampling along the spatial direction /?. The problem to be solved is thus one of inverting the 
spectral aliasing (i.e. the over-resolution) using a finite number of discrete data provided by a complex 
system. The solution proposed here is based on precise modelling of the instrument and, in particular, 
the integral equations containing the continuous variables (a, /5 and A) of the optics and sensing system. 
The model input is naturally a function of these continuous variables (\>(a^ P, A) and the output is a finite 
set y of discrete data items. 

'in this paper, the spatial dimensions are angles in radian 
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The approach used for solving the inverse problem, i.e. reconstructing an object having three continuous 
variables from the discrete data 

• comes within the framework of regularization by penalization; 

• uses a semi-parametric format where the object is decomposed into a family of functions. 

There is a multitude of families of functions available (possibly forming a basis of the chosen functional 
space). The most noteworthy are Shannon, Fourier, wavelet and pixel-indicator families or those of spline, 
Gaussian Kaiser-Bessel, etc. Work on 3D tomographic reconstruction has used a family of Kaiser-Bessel 
functions having spherical symmetry in order to calculate the projections more efficiently 121, 131, El, Q. 
In a different domain, the signal processing community has been working on the reconstruction of over- 
resolved images from a series of low resolution images ISJ. A generic direct model can be described |6|, 
starting with a continuous scene, to which are applied k shift or deformation operators including at least 
one translation. This step gives k deformed, high-resolution images. A convolution operator modelling 
the optics and sensor cells is then applied to each of the images. After subsampling, the k low-resolution 
images that constitute the data are obtained. Recent work on the direct model has mainly concerned 
modelling the shift by introducing a rotation of the image Q, lH and a magnifying factor lOj. Other 
works have modelled the shift during sensor integration by modifying the convolution operator ITOl . To 
the best of our knowledge, in most works, the initial discretization step is performed on pixel indicators 
EH, lEl, HOl, HI, im. on this point, a noteworthy contribution has been made by P. Vandewalle 
et al. who discretize the scene on a truncated discrete Fourier basis HH. However, their decomposition 
tends to make the images periodic leading to create artefacts on the image side. Thus we have decided 
not to use this approach. Recently, the problem of X-ray imaging spectroscopy has been solved in the 
Fourier space |14], but each spectral component has been estimated independantly. 

The two major contributions of our paper are (1) the modelling of the measurement system as a whole 
with continuous variables and (2) the continuous variable decomposition of the three dimensional object 
over a family of Gaussian functions. Modelling with continuous variables enables a faithful description 
to be made of the physical phenomena involved in the acquisition and avoids to carry out any prior data 
interpolation. In our case, computing the model output requires six integrals (two for the response of the 
optics, two for the grating response, and two for the sensor integration) and the choice of a Gaussian 
family allows five of these six integrals to be explicitly stated. Our paper is organised as follows: 

Part In] describes the continuous model of the instrument comprising: the diffraction at the aperture, the 
truncation by the slit, the response of the grating, the distortion of the light flux, the sensor integration 
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and the scanning of the sky. In part Hill the object with continuous variables is decomposed over a family 
of Gaussian functions. The aperture and grating responses are approximated by Gaussian functions. This 
part concludes with the obtention of a precise, efficient model of the measuring system. The inverse 
problem is solved in a regularized framework in part |IVl Finally, part |V] gives an evaluation of the 
resolving power of the method and a comparison with a standard data co-addition method using both 
simulated and real data. 

II. Continuous direct model 

The aim of the instrument model is to reproduce the data, y, acquired by the spectral imager from 
a flux 4>{a, (3, A) of incoherent light. Fig. [T] illustrates the instrument model for one acquisition (the 
telescope remains stationary). To simplify, we present the scanning procedure in section III-DI First, we 
have the response of the primary mirror (aperture), which corresponds to a convolution. Second, there is 
a truncation due to a rectugular slit. Third, a grating disperses the light. Finnally, the sensor integration 
provides the discrete data y. Distortion of the luminous flux is modelled in the sensor integration. 



Sky 



Aperture 
diifraction 



Distortion 



Sensor integration 



4>f 



Slit response 



Diffraction grating 



Model out- 
put y{i,j) 



Fig. I. Block diagram of the direct model for one acquisition: from a continuously defined sky (/> to a discrete output y describing 
the data. The flux is a convolution of the flux and the PSF of the primary miror. <l>f is truncated by a rectangular slit and 
is dispersed by the gratting. Finally, the sensor provide a discrete output y. 



A. Aperture diffraction 

Under some hypotheses, the propagation of a light wave which passes through an aperture is determined 
by Fresnel diffraction i fTSl and the result in the focal plane is a convolution of the input flux (f) with 
the Point Spread Function (PSF) ha illustrated in Fig. |2] for a circular aperture. This PSF, which is a low 
pass filter, has a width proportional to the wavelength of the incident flux. For a circular aperture, it can 
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be written: 



ha{a,(3,X) = A 



Ji(7rZ?v'«' + /3VA) 



(1) 



where Ji is the first order Bessel function of the first kind, A is an amplitude factor and D is the diameter 
of the mirror. 




A] = 15/(rn 
A2 = 7/im 



i 22'^ 



D 



Fig. 2. Profile of an Airy disk (PSF for a circular aperture) for two wavelengths. 



The flux in the focal plane, (pj, is written in integral form: 

/?', A) = [[ (l){a,p,X)haia-a',P-P',X)dadp 

J Ja,/3 



(2) 



B. Slit and diffraction grating 

a) Slit: Ideally, the slit and grating enable the dispersion of the wavelengths in spatial dimension (3 
previously "suppressed" by the slit (see Fig. [3]). In practice, the slit cannot be infinitely narrow because 
the flux would be zero. The slit thus has a width 7 of about two pixels. 

b) Diffraction grating: Ideally, the grating gives a diffracted wave with an output angle 9 linearly 
dependent on the wavelength A (see Fig. [3]). In a more accurate model, the dependencies become more 
complex. Let us introduce a variable u in order to define an invariant response of the system lfT6ll . 

sin 6 — sin (3' sin — (3' 
u = (3) 

where (3' is the angle of incidence of the wave on the grating , and < 7/2 where 7 is the angular 
slit width (5.6 arcseconds). The response of the grating centred on mode m (m = 0, 1, . . . ) can, with 
some approximations, be written as the square of a cardinal sine centred on m/a [,16] . 

hr{e,P',X) = Bsmc'^{TTL{u-m/a)) (4) 
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Grating 




Detector 



Fig. 3. Optical scheme of IRS instrument: the sht is on the focal plan of the telescope Spitzer. The gratting disperses the 
light and the detector collectes the dispersed flux 



where L is the width of the grating and a the grid step (distance between two grooves). This response 
centred on the first mode (m = 1) is plotted in Fig. |4l 




1/a 



Fig. 4. Diffraction grating response. The grey curve corresponds to the first mode. In reality, the response width is smaller. 

As the flux is an incoherent light source, the expression for the signal at the output of the grating is 
written in the form of an integral over [3' and A: 



^r{a',e)= / (l)f{a',f3',X)hr{P',X,9)d/3'dX 

J\ ^l/3'|<7/2 



(5) 



where / is the slit width. 
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C. Sensor integration 

Once the flux has passed through the grating and the wavelengths have been dispersed according to 
6, the Ught flux is focused on the sensor composed of square detectors. The sensor is simply modelled 
by integrating the flux (j)^ on square areas of side d. The flux is integrated along the direction a, which 
is not modified by the diffraction grating, and the dimension 6, a combination of (3 and A, to obtain the 
discrete values 

[■(i+l)d [■(j+l)d+e1. 

y{i,j) = / / Mc^',0)da'd9. (6) 

J id J jd+ej^ 

The integration limits are modified by the terms in order to take into account the data distortion 
as illustrated in Fig. [5] 



e 

Fig. 5. Modelling the distortion: sensor integration limits are shifted according to the dimensions a . 

D. Scanning procedure of the sky 

In a direction parallel to the slit width a scanning procedure (illustrated Fig. O is applied. This scanning 
procedure is composed of Q acquisitions. Between the first and the q*^ acquisitions, the instrument is 
moved by ^a{(l) (resp. A^(g)) in the direction a (resp. (3). To taking into account the motion of the 
instrument, we substitue 0(a, /?, A) for (/)(a — /S.a{q) , (3 — /S.p{q) , A) in the previous equations. In practice, 
we fix the a axis in the direction of the slit and the j3 axis perpendicular to the slit (see Fig. |6ll. In 
consequence, Aq,((7) is equal to zero. 

E. Complete model 

By combining expression ([T|), (O, dS),® and Q, we obtain a continuous direct model in the form 
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a J/3 

where ^ is a scale factor. 

The equation dTJ can rewritten: 



/ 

id Jjd+e]. J X J\l3'\<'r/2 

a - Aa{q), P - A/siq), \)ha{a - a' , (3 - /?', A)dad/? 



hr{e,(3',X)dp'dXda'd9 (7) 



y{^,j,q)= / 4>{a,P,X)h't^f{a,P,X)dadpd\ (8) 

J a J/3 J\ 

with 

hl'^f{a,P,X)=A / / 

J id Jjd+elj J-lfi 

ha{a -a'- A,(g), /3 - /3' - A^(g), X)hr{9, /3' , X)da'd9dP' (9) 

We have been developed a model relying the continuous sky /3, A) and discrete data y. Our 
model is linear not-shift-invariant, because the aperture response and the grating response depend on 
the wavelength. 

III. Decomposition over a family and Gaussian approximation 

In the previous part, we have seen that obtaining the output from the model requires the six integrals 
of equation (|7]) to be calculated. The estimation of in L^(M) by inversion of this model is quite tricky, 
so we prefer to decompose the object over a family of functions. As we can see in the introduction, a lot 
of such decomposition functions can be used. The most traditional are Fourier bases, wavelets, cardinal 
sines, splines and pixel indicators. The choice does not have any great influence on the final result if 
the continuous object is decomposed over a sufficiently large number of functions. We therefore chose 
our decomposition functions in such a way as to reduce the computing time for the instrument model. 
First, we chose the a axis in the direction of the slit and the /3 axis perpendicular to the slit (see Fig. [6]). 
Second, we have two spatial variables (a, /3) and one spectral variable A, so to simplify the calculus, we 
chose decomposition functions that are separable into (a, /?) and A. Thrid, the object is convolved by the 
response of the optics, which has circular symmetry. So we choose functions possessing the same circular 
symmetry in order to make this calculation explicit. Finally, the slit and the grating have an impact in 
the P direction only ([5]), which motivates us to choose functions that are separable into a and /3. These 
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considerations led us to choose Gaussian functions along the spatial directions. Finally the complexity 
of the A dependence encouraged us to choose Dirac impulses for the spectral direction. 

A. Decomposition over a family of Gaussian functions 

The flux is a continous function decomposed over a family of separable functions: 

k I p 

U{a - kT^) - ITp) r(A - pTx) 

k I p 

where x{k,l,p) are the decomposition coefficients, Ta, Tp and Tx are the sampling steps, and with: 



r(A) = S{\) (12) 

With such decomposition, the inverse problem becomes one of estimating a finite number of coefficients 
x{k,l,p) from discrete data y{i,j,q). By combining equations ([8]l and ([TOl l. we obtain: 



yihj,Q)='^'^'^x{k,l,p) / j '^k,i,p{a,l3,X) 



k I p 



a J/3 JX 



hl'^i'^{a,(3,X)dadpdX (13) 

If the y{i,j,q) and x{k,l,p) are gathered in vectors y and ao respectively, the equation ([T3l ) can be 
formalized as a vector matrix product. 

y = Hx (14) 

with each component of the matrix H is calculated using the integral part of the equation ([T3] ). The 
n{k, /,^')-th column of matrix H constitutes the output when the model is used with the n-th decompo- 
sition function i"^k,i,p)- The model output for '^k,i,p is calculated in the next two sections. 

^In this paper, we use the following convention: bold, lower-case variables represent vectors and bold, upper-case variables 
represent matrices. 
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B. Impulse responses approximated by Gaussian functions 

1 ) Approximation of the PSF : Equation Q comes down to convolutions of a squared Bessel function 
and Gaussians. This integral is not explicit and, in order to carry out the calculations, the PSF is 
approximated by a Gaussian 



with a standard deviation ax depending on the wavelength. Indeed, the Bessel functions cross zero at 
the first time in 1.22X/D. ax is determined numerically by minimizing the quadratic error between the 
Gaussian kernel and the squared Bessel function, which gives for our instrument ax ~ A/2. The relative 
quadratic error errL2 = ||Bessel — Gaussian^ / ||Bessel||2 is equal to 0.15% for our instrument. If we 
caculate the relative absolute error err^i = ||Bessel — Gaussian||i / ||Bessel||i, we obtain 5 %. We can 
conclude that most of the energy of the squared bessel function is localized in the primary lobe. Another 
advantage of using the Gaussian approximation is that the convolution kernel is separable into a and /?. 
Finally, the result of the convolution of two Gaussian functions is a standard one and is also a Gaussian: 



2 ) Approximation of the grating response: The presence of the slit means that integral dSjl is bounded 
over /?' and is not easily calculable. Since the preceding expressions use Gaussian functions, we approx- 
imate the squared cardinal sine by a Gaussian to make the calculations easier: 



as is determined numerically by minimizing the quadratic error between the Gaussian kernel and the 
squared cardinal sine, which gives for our instrument as ~ 25.5 m~^. The relative errors made are larger 
than the bessel case = 0.43%,errj^i = 10.7%), but this Gaussian approximation of the grating 

response allows the flux (pr coming out of the grating to be known explicitly. 

The error introduced here is larger than for the Gaussian approximation of the PSF described in the 
previous section. However, our goal is to have a good model of the spatial dimension of the array. 
Furthermore, with respect to the current method, the fact of taking the response of the grating into 
consideration, even as an approximation, is already a strong improvement. 

February 11, 2009 DRAFT 





(15) 





(16) 



SUBMITTED TO IEEE JOURNAL OF SELECTED TOPICS IN SIGNAL PROCESSING 



11 



Ma',d) = I F'^ I [ n{a - - Pi - Ap{q)) 

r(A - Xp)ha{a -a, 15- /?', \)hr{e, f3' , A) 
dad/5d/3'dA 



with 



, (a' - akf 
^exp — -T-n ^ I exp 



sin( 



A 

V 

erf(a) 



erf 



erf 



2S2 

7/2 -/u 



4vry (a2 + a2)(a2 + a2 + A2a2) 
^ + A + A/3((z) 



smi 



(a2+a2) + A2a2 



(17) 



In equation (fTTl) . it can be seen that c/i^ is separable into a' and ^. Let us introduce the functions / and 
g such that: 

Ma',e) = Af{a')g{9) (18) 



C. Sensor integration 
First, we calculate the sensor integration in the a' direction. 

r(j+l)d+e% 



/(a')da' 



erf 



{j + l)d + e| - a k 



erf 



^2K + 



(19) 



with /C = J^(cjf + cj2)/2 
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The integral of g is calculated numerically as the presence of erf functions in equation (fTTl) does not 
allow analytical calculations. 

We obtain the expression for the n-th column of matrix H, which now contains only a single integral: 

f{j+l)c( 



y{i,j,q)=AIC g{9)d9 

J id 



U + l)d + e?- -ak\ J jd + ej, - Uk 



(20) 



erf I ^ - erf ' 

Using expression (l20l ). the elements of matrix H are pre-computed relatively rapidly. Thanks to the 
sparsity of the matrix H to calculate the model output of Eq. ([141) . 

IV. Inversion 

The previous sections build the relationship (fT4l ) between the object coefficients and the data: it 
describes a complex instrumental model but remains linear. The problem of input (sky) reconstruction is 
a typical inverse problem and the literarure on the suject is abundant. 

The proposed inversion method resorts to linear processing. It is based on conventional approaches 
described in books such fTTl . iflSl or more recently [19|. In this framework, the reader may also 
consider ||20l . 11211 for inversion based on specific decomposition. These methods rely on a quadratic 
criterion 

J{x) = \\y - Hx\\^ + ljLajj\\Dapx\\^ + ^,x\\Dxx\\^ . (21) 

It involves a least squares term and two penality terms concerning the differences between neighbouring 
coefficients: one for the two spatial dimensions and one for the spectral dimension. They are weighted by 
^ap and /Ua, respectively. The estimate x is chosen as the minimizer of this criterion. It is thus explicit 
and linear with respect to the data: 

x={H'H + fi^pDl^^D^f, + fixD\Dx) H'y (22) 

and depends on the two regularization parameters and jj-x- 



Remark 1 — This estimator can be interpreted in a Bayesian framework l\22\l based on Gaussian 
models for the errors and the object. As far as the errors are concerned, the model is a white noise. 
As far as the object is concerned, the model is correlated and the inverse of the correlation matrix is 
proportional to fj,ai3DapDaf3 + fJ-xD^Dx, i.e., it is a Gauss markov field. In this framework, the estimate 
maximizes the a posteriori law. 
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Remark 2 — Many works in the field of over-resolved reconstruction concern edge preserving pri- 
ors / l23l/ . Aiil/ . /[7|/, A24I/ . [1], / f72l/ . /?i owr application here, smooth interstellar dust clouds are under 
study, so preservation of edges is not appropriate. For the sake of simplicity of implementation, we chose 
a Gaussian object prior. 

The minimizer x given by relation (l22l) is explicit but, in practice, it cannot be calculated on standard 
computers, because the matrix to be inverted is too large. The solution x is therefore computed by 
a numerical optimization algorithm. Practically, the optimization relies on a standard gradient descent 
algorithm |[25l . ll26ll . More precisely, the direction descent is a approximate conjugate gradient direction 
|[27l and the optimal step of descent is used. Finally, we initialise the method with zero (x = 0). 

V. Results 

As we have presented in part |lll] the a and j3 axis are fix (see Fig. [6l right). The real data is composed 
of 23 acquisitions having a spatial dimension a' and a spectral dimension of wavelength between 7.4 
and 15.3 jini (each acquisition is an image composed of 38 x 128 detector cells, see Fig. [6l left). Between 
two acquisitions, the instrument is moved by half a slit width in the /3 direction. Fig. [6l right, shows the 
scanning procedure applied to the Horsehead nebula |[28l . 



15,3 /im 




Fig. 6. Acquisition. Tlie left-hand image represents tlie data acquired for one pointing position: the vertical axis shows the 
spectral dimension 6 and the horizontal axis is the spatial dimension a' . The slit is represented schematically by a rectangle in 
the middle. The right-hand image illustrates the scanning strategy in the j3 direction. 



Our results (Fig. |8(b)[ solid line on Fig. |9] and Fig. |10(b)[ ) can be compared with those obtained with the 
conventional processing (Fig. 8(c) dotted line on Fig. |9] and Fig. |10(a)[ ). For the conventional processing 
(described in Compiegne et al. 2007 l28l) an image of the slit is simply extracted for each wavelength 
from the data taken after each acquisition (e.g. left panel of Fig. O and projected and co-added on the 
output sky image, without any description of the instrument properties. 
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A. Simulated data 

In our first experiment, we reconstruct data simulated using our direct model. We choose an object 
with the same spatial morphology and the same spectral content as the Horsehead nebula (see Fig. |8(a)[ ). 
However, in order to tune the regularization coefficient, we perform a large number of reconstructions. 
Thus, we need to simulate a problem smaller than in our real case. The data are composed of 14 
acquisitions, and the virtual detector contained 18 x 40 pixels. We choose to reconstruct a volume with 
15870 gaussians distributed on a cartesian grid 23 x 23 x 30. Finally, we add to the output of the model 
a white Gaussian noise with the same variance as the real data. 

The results contain a set of 30 images (see Fig. [T]). Fig. |8(b)| and solid line on Fig|9] illustrate our result 
for one wavelength (8.27 fim) and one pixel, respectively. The image computed with our method (Fig. 
|8(b)| ) appears comparable to the true image (Fig. |8(a)| ), while the image computed with the conventional 
processing (Fig. |8(c)] l is smoother. A comparison of solid line and dotted line in Fig. |9] clearly also shows 
that our method provides a spectrum comparable to the true spectrum, while the peaks obtained with the 
conventional processing are too broad. 

Our sky estimation depends on the regularization coefficients ^1^/3 and fix- We tune this parameters 
by minimizing numerically the quadratic error between the estimated object and the real object were 
selected. In this experiment we obtain fia/s = 0.01, fix = 0.005. 

B. Real data 

Eeal data contain 23 acquisitions composed of 38 x 128 values. To obtain a over-resolved reconstruction, 
we describe our volume with 587264 gaussians destributes on a cartesian grid 74 x 62 x 128. The spatial 
(a, (3) sampling step is equal to a quater slit width, and the spectral dimension is uniformly sampled 
between the wavelength 7.4 and 15.3 fim. The reconstruction is computed after setting the regularization 
coefficients fi^/s and fix empirically. Too low a value for these coefficients produces an unstable method 
and a quasi explosive reconstruction. Too high a value produces images that are visibly too smooth. A 
compromise found by trial and error led us to fi^js = 0.3 and fix = 0.7. The ratio between fi^is and fix 
is also based on our simulation. However, we cannot compare the regularization coefficients between the 
simulated and the real case, since the size of the problem modifies the weigth of the norm in the Eq. 
(|2T]) . Pratically, we take large value for the regularization coefficients, and we gradually reduce the value 
up that we are seeing noise. 

Our results (Fig. |10(b)| and |ll(b)[ ) can be compared with those obtained with (Fig. |10(a)| and |ll(a)[ 
from [28]). A comparison of Fig. 10(a) and 1 10(b)] clearly shows that our approach provides more resolved 
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Fig. 7. Set of 30 images of our reconstruction from simulated data. Each image corresponds to one wavelength for 7.4 to 
9.2 fim with a step of 0.062 fim. 




(a) (b) (c) 



Fig. 8. Image at A = 8.27 /im:, (a) simulated sky, (b) image estimated by our method, (c) image estimated by a conventional 
method. 
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Fig. 9. Spectrum of a pixel. The curves abscissa is the wavelenght in meters: solid line: simulated sky, dashed line: our method, 
dotted line: conventional method. 
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images that bring out more structures than the conventional approach. Note, in particular, the separation of 
the two filaments on the left part of the Fig. |10(b)| obtained with our method, which remains invisible after 
conventional processing. For comparison. Fig. |10(c)| shows the same object observed with the Infrared 
Array Camera (IRAC) of the Spitzer Space Telescope which has a better native resolution since it observes 
at a shorter wavelength (4.5 /im). Here the same structures are observed, providing a strong argument in 
favour of the reality of the results provided by our method. 

A more precise analysis is done in section IV-CI It provides a quantitative evaluation of the resolution. 

Finally, the spectra reconstructed by our method (Fig. |1 l(b)[ ) have a resolution slightly better than the 



one reconstructed by the conventional method (Fig. |1 l(a)| ). The peaks characterizing the observed matter 
(gas and dust) are well positioned, narrower and with a greater amplitude. However, ringing effects appear 
at the bases of the strongest peaks (Fig. |ll(b)[ ). They could be explained by an over-evaluation of the 
width (Ts of the response of the grating, or by the gaussian approximation. 





(a) 



(b) 




(c) 



Fig. 10. Reconstruction of a slcy representing the Horsehead nebula: (a) image estimated at 11.37 fj,m by the conventional 
method (b) image estimated at 11.37 by our method, (c) image obtained with the Infrared Array Camera IRAC on board 
the Spitzer Space Telescope having better resolution at 4.5 iim.. 




7.4 8.6 9.9 11.1 12.4 13.6 14.8 (1 m 7.4 8.6 9.9 11.1 12.4 13.6 14.8 n m 

(a) (b) 

Fig. 11. Spectrum of the center point of the Fig. |10(b")l (a) spectrum estimated with the conventional method (b) spectrum 
estimated with our method. 
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C. Study of resolving power of our approach 

This part is devoted to numerical quantification of the gain in angular resolution provided with 
our method, using the Rayleigh criterion, which is frequently used by astrophysicists: for the smaller 
resolvable detail, the first minimum of the image of one point source coincides with the maximum of 
another In practice, two point sources with the same intensity and a flat strectrum are considered to 
be separated if the minimal flux between the two peaks is lower than 0.9 times the flux at the peak 
positions. The resolution is studied in the (5 direction only as this is the direction in which the subslit 
scan is performed. 

Two point sources are injected, at positions f3i and (32, respectively (see Fig. [131 top). The corresponding 
data are simulated, and the reconstruction (/>(/?) is performed. As explained above, the two point sources 
are considered to be separated if (j){[l5i + (52] /2) < 0.9x The resolution is defined as the difference 

(5 = /32 — /?i at which the two point sources start to be separated. 

Point sources are simulated for a set of differences 5 between 2.4 and 5.4 arcseconds and simulations 
are performed in the configuration of the real data (signal to noise ratio, energy of the data). Moreover, 
we use the regularization parameters and determined in section IV-AI A number of reconstructions 
has been obtained. The ratio between the values of the reconstructed function at /?i and {(5i + /32)/2 is 
calculated as a function of the difference 5 between the two peaks. Results are shown in Fig. [T2l 

The computed resolution is 3.4 arcseconds (see Fig. H^a)) and 5 arcseconds (see Fig. fl^ b)) for our 
method and the conventional method, respectively. Fig. [13] illustrates this gain in angular resolution. In 
the left column on Fig. [13] {5 = 3.4) corresponds to the limit of resolution of our method. In this case, 
the peak is not separated with the conventional method (Fig. [13] (d)). In the middle column on Fig. [13] 
our algorithm clearly separates the peak (Fig. [T3lh)) and not the conventional method (Fig. [T3le)). In the 
right column, we observe a stain with our method is smaller than the conventional method. Our method 
increases the resolution by a factor 1.5. 

VI. Conclusions 

We have developed an original method for reconstructing the over-resolved 3D sky from data provided 
by the IRS instrument. This method is based on: 

1) a continuous variable model of the instrument based on a precise integral physical description, 

2) a decomposition of the continuous variable object over a family of Gaussian functions, which results 
in a linear, semi-parametric relationship. 
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(b) 

Fig. 12. Resolution of our method: the curve represents the ratio of the intensity at one peak to the intensity between the two 
peaks as a function of the distance between the peaks in arcseconds. The resolution is read at crossing of this curve and the 
dotted line (the ratio is 0.9). (a) Results obtained with our method, (b) Results obtained with the conventional method 




3) an inversion in the framework of deterministic regularization based on a quadratic criterion mini- 
mized by a gradient algorithm. 

The first results on real data show that we are able to evidence spatial structures not detectable using 
conventional methods. The spatial resolution is improved by a factor 1.5. This factor should increase 
using data with a motion between two acquisitions smaller than the half a sUt width. 

In the future, we plan to design highly efficient processing tools using our approach in particular for 
the systematic processing of the data which will be taken with the next generation of infrared to milimeter 
space observatory (Herschel, Planck, ...). 
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